///*
// * Measure.cpp
// *
// *  Created on: Apr 20, 2010
// *      Author: snirgaz
// */
//
#include "Measure.h"
//#include <array>
//using namespace std;
//
//Measure::Measure() {
//	closedMeasure_ = openMeasure_ = false;
//	numOfClosed_ = numOfOpen_ = 0;
//	VectorObs<double>::InitData initDataVec = { QmcParams::numOfBins,
//			QmcParams::M, QmcParams::binSize };
//	ScalarObs<double>::InitData initDataScalar = { QmcParams::numOfBins,
//			QmcParams::binSize };
//	// Generate FFT Buffer
//	if (QmcParams::toMeasure.Density || QmcParams::toMeasure.JqJq
//			|| QmcParams::toMeasure.DensitySym) {
//		array<int, 1> tauSize = { QmcParams::M };
//		fftPlan_.init(&sfSum_[0], &sfSum_[0], 1, &(tauSize[0]));
//	}
//	if (QmcParams::toMeasure.Density || QmcParams::toMeasure.DensitySym) {
//		closedMeasure_ = true;
//		sfData_.resize(QmcParams::NumP * QmcParams::M);
//		sfSum_.resize(QmcParams::M);
//		//sfData_ = ippsMalloc_64fc(QmcParams::NumP * QmcParams::M);
//		//sfSum_ = ippsMalloc_64fc(QmcParams::M);
//
//	}
//	// Structure Factor Measurements
//	if (QmcParams::toMeasure.Density || QmcParams::toMeasure.DensitySym) {
//		array<int, SIM_DIM> sizes;
//		sizes.fill(QmcParams::numJqJqPoints);
//		sfqw_.init(sizes, initDataVec);
//		sf_.init(initDataScalar);
//		sfPow_.resize(QmcParams::M);
//		std::fill(sfPow_.begin(), sfPow_.end(), 0.);
//	}
//	if (QmcParams::toMeasure.DensitySym) {
//		array<int, 1> sizes;
//		sizes[0] = QmcParams::numOfKvecSym;
//		VectorObs<double>::InitData initDataVec = { QmcParams::numOfBins,
//				QmcParams::M / 2 + 1, QmcParams::binSize };
//		sfqwSym_.init(sizes, initDataVec);
//		kVec_.resize(QmcParams::numOfKvecSym);
//		int G = int(QmcParams::period), diff = -QmcParams::numOfKvecSym / 2 + 1;
//		for (auto itk = kVec_.begin(); itk != kVec_.end(); itk++) {
//			array<int, SIM_DIM> pos;
//			for (int d = 0; d < SIM_DIM; d++) {
//				pos[d] = 0;
//				pos[(d + 1) % SIM_DIM] = G + diff;
//				pos[(d + 2) % SIM_DIM] = G + diff;
//				itk->push_back(pos);
//				pos[(d + 2) % SIM_DIM] = -(G + diff);
//				itk->push_back(pos);
//			}
//			diff++;
//		}
//		closedMeasure_ = true;
//	}
//#if SIM_DIM==2
//	if (QmcParams::toMeasure.DensityPhonon) {
//		int kmax = int(QmcParams::period);
//		int size = kmax * kmax * QmcParams::M;
//		double dx = 1. / double(kmax);
//		MKL_LONG status = 0;
//		MKL_LONG N[3];
//		N[0] = kmax;
//		N[1] = kmax;
//		N[2] = QmcParams::M;
//		closedMeasure_ = true;
//		boost::array<int, 3> sizes;
//		sizes[0] = sizes[1] = kmax;
//		sizes[2] = QmcParams::M / 2 + 1;
//		sfPhononqwObs_.resize(2);
//		for (int d = 0; d < QmcParams::dim; d++) {
//			deltau_.push_back(new Array3D<double>(kmax, kmax, QmcParams::M));
//			deltauqw_.push_back(
//					new Array3D<std::complex<double>>(kmax, kmax,
//							QmcParams::M / 2 + 1));
//			//	sfPhononqw_.push_back(
//			//		new Array3D<double>(kmax, kmax, QmcParams::M / 2 + 1));
//			Plans_.push_back(
//					new fftPlan(deltau_[d].data(), deltauqw_[d].data(), 3, N));
//			sfPhononqwObs_[d].init(sizes, initDataScalar);
//		}
//
//	}
//#endif
//	if (QmcParams::toMeasure.Angle) {
//		closedMeasure_ = true;
//		int size = 1;
//		for (int d = 0; d < QmcParams::dim; d++)
//			size *= QmcParams::numJqJqPoints;
//		array<int, SIM_DIM> sizes;
//		sizes.fill(QmcParams::numJqJqPoints);
//		angleqw_.init(sizes, initDataVec);
//		angleZeroQ_.init(initDataScalar);
//		nnIndices_.resize(
//				boost::extents[QmcParams::M][QmcParams::NumP][QmcParams::NumP]);
//		angle_.resize(QmcParams::NumP * QmcParams::M);
//		//angle_ = ippsMalloc_64fc(QmcParams::NumP * QmcParams::M);
//	}
//	// Green's Function (Density Matrix rho(0,r)
//	if (QmcParams::toMeasure.Rho) {
//		openMeasure_ = true;
//		array<int, SIM_DIM> sizes;
//		std::fill(sizes.begin(), sizes.end(), QmcParams::numJqJqPoints);
//		rhoQ_.init(sizes, initDataScalar);
//		// Momentum space
//		//rhoQ_ = ippsMalloc_64f(size);
//		//ippsZero_64f(rhoQ_, size);
//	}
//	if (QmcParams::toMeasure.Winding || QmcParams::toMeasure.JqJq) {
//		velocity_.resize(QmcParams::NumP, QmcParams::M, SIM_DIM);
//	}
//	if (QmcParams::toMeasure.DensitySym || QmcParams::toMeasure.Density
//			|| QmcParams::toMeasure.JqJq) {
//		expData_.resize(QmcParams::NumP * QmcParams::M);
//	}
//	if (QmcParams::toMeasure.Winding) {
//		closedMeasure_ = true;
//		windSqr_.init(initDataScalar);
////		int length = 2 * QmcParams::NumP + 1;
////		windDist_.resize(length);
////		std::fill(windDist_.begin(), windDist_.end(), 0);
//	}
////	if (QmcParams::toMeasure.JqJq) {
////		closedMeasure_ = true;
////		int numOfBeads = WL::size();
////		boost::array<int, SIM_DIM> sizes;
////		std::fill(sizes.begin(), sizes.end(), QmcParams::numJqJqPoints);
////		JqJq_.init(sizes, QmcParams::numOfBins, QmcParams::M);
////		averagePos_ = new WorldLinesData(QmcParams::NumP, QmcParams::M,
////				SIM_DIM);
////		velocityComplex_ = ippsMalloc_64fc(numOfBeads);
////		Jq_ = ippsMalloc_64fc(numOfBeads);
////		JqSum_ = ippsMalloc_64fc(QmcParams::M * SIM_DIM);
////		JqPow_ = ippsMalloc_64f(QmcParams::M);
////	}
//	if (QmcParams::toMeasure.Energy) {
//		EK_.init(initDataScalar);
//		EPV_.init(initDataScalar);
//		EPF_.init(initDataScalar);
//		E_.init(initDataScalar);
//		closedMeasure_ = true;
//	}
////	if (QmcParams::toMeasure.Hexatic) {
////dr_	=
////}
//	numOfClosed_ = 0;
//	numOfOpen_ = 0;
//	//	ippGetCpuFreqMhz(&mhz_);
//}
//Measure::~Measure() {
//	if (QmcParams::toMeasure.Density || QmcParams::toMeasure.JqJq
//			|| QmcParams::toMeasure.DensitySym) {
//		//	ippsFree(fftBuffer_);
//		//	ippsDFTFree_C_64fc(ppDFTSpec_);
//		//	ippsFree(expData_);
//	}
//	if (QmcParams::toMeasure.Density || QmcParams::toMeasure.DensitySym) {
//		//	ippsFree(sfData_);
//		//	ippsFree(sfSum_);
//	}
//	if (QmcParams::toMeasure.Angle) {
//		//	ippsFree(angle_);
//	}
//	if (QmcParams::toMeasure.Rho) {
//
//	}
//	if (QmcParams::toMeasure.Winding || QmcParams::toMeasure.JqJq) {
//		//	delete velocity_;
//	}
//	if (QmcParams::toMeasure.JqJq) {
////		ippsFree(JqJq_);
////		ippsFree(JqSum_);
////		ippsFree(JqPow_);
////		ippsFree(Jq_);
////		delete averagePos_;
////		ippsFree(velocityComplex_);
//	}
//}
//void Measure::measureDensity() {
//	int kmax = QmcParams::numJqJqPoints;
//	vector<complex<double>>::iterator sfDataPos, sfSumPos;
//	int dimSize = WL::numOfBeads();
//	ConfDate& wlData = WL::data();
//	DoubleVec sf;
//	DoublePointerVec wlDim;
//	double norm = 1
//			/ double(
//					QmcParams::NumP * QmcParams::NumP * QmcParams::M
//							* QmcParams::M);
//	double normet = 1
//			/ double(QmcParams::NumP * QmcParams::NumP * QmcParams::M);
//	WL::getDimP(wlDim);
//	for (auto it = ++sfqw_.begin(); it != sfqw_.end(); it++) {
//		array<int, SIM_DIM> kVec = it.getIndex();
//#if SIM_DIM==2
//		expData_[0:dimSize]=2*PI*(
//				double(kVec[0])*wlDim[0][0:dimSize]
//				+double(kVec[1])*wlDim[1][0:dimSize]);
//
//#endif
//#if SIM_DIM==3
//		for (int i=0; i< dimSize;i++) {
//			sfData_[i]=polar(1.,2*PI*(
//							double(kVec[0])*wlDim[0][i]
//							+double(kVec[1])*wlDim[1][i]
//							+double(kVec[2])*wlDim[2][i]));
//		}
//#endif
//		// ippsCIS_64fc_A26(expData_, sfData_, dimSize);
//		//sfDataPos = sfData_;
//		//sfSumPos = sfSum_;
//		sfDataPos = sfData_.begin();
//		sfSumPos = sfSum_.begin();
//		double sfSlice = 0;
//		for (int m = 0; m < QmcParams::M; m++) {
//			//ippsSum_64fc(sfDataPos, QmcParams::NumP, sfSumPos);
//			(*sfSumPos) = accumulate(sfDataPos, sfDataPos + QmcParams::NumP,
//					complex<double>(0., 0.));
//			sfSlice += std::norm(*sfSumPos);
//			//(pow((*sfSumPos).re, 2) + pow((*sfSumPos).im, 2));
//			sfDataPos += QmcParams::NumP;
//			sfSumPos++;
//		}
//		fftPlan_.execute();
//		//ippsDFTFwd_CToC_64fc(sfSum_, sfSum_, ppDFTSpec_, fftBuffer_);
//		sf.push_back(sfSlice * normet);
//		// Power Spectrum
//		transform(sfSum_.begin(), sfSum_.end(), sfPow_.begin(),
//				[](complex<double> sf)->double {return std::norm(sf);});
//		// ippsPowerSpectr_64fc(sfSum_, &(sfPow_.data()[0]), QmcParams::M);
//		sfPow_ *= norm;
//		(*it) << (sfPow_);
//	}
//	sf_ << (*std::max_element(sf.begin(), sf.end()));
//}
//void Measure::measureDensitySym() {
//	vector<complex<double>>::iterator sfDataPos, sfSumPos;
//	int dimSize = WL::numOfBeads();
//	DoublePointerVec wlDim;
//	double norm = 1.
//			/ double(
//					QmcParams::NumP * QmcParams::NumP * QmcParams::M
//							* QmcParams::M);
//	WL::getDimP(wlDim);
//	auto itObs = sfqwSym_.begin();
//	for (auto itKVec = kVec_.begin(); itKVec != kVec_.end(); itKVec++) {
//		double AvNumOfKVev = 1. / int(itKVec->size());
//		sfPow_.resize(QmcParams::M / 2 + 1);
//		std::fill(sfPow_.begin(), sfPow_.end(), 0.);
//		for (auto itKVecC = itKVec->begin(); itKVecC != itKVec->end();
//				itKVecC++) {
//			array<int, SIM_DIM> kVec = *itKVecC;
//#if SIM_DIM==2
//			expData_[0:dimSize]=2*PI*(
//					double(kVec[0])*wlDim[0][0:dimSize]
//					+double(kVec[1])*wlDim[1][0:dimSize]);
//
//#endif
//#if SIM_DIM==3
//			for (int i=0; i< dimSize;i++) {
//				sfData_[i]=polar(1.,2*PI*(
//								double(kVec[0])*wlDim[0][i]
//								+double(kVec[1])*wlDim[1][i]
//								+double(kVec[2])*wlDim[2][i]));
//			}
//#endif
//			//ippsCIS_64fc_A26(expData_, sfData_, dimSize);
//			sfDataPos = sfData_.begin();
//			sfSumPos = sfSum_.begin();
//			double sfSlice = 0;
//			for (int m = 0; m < QmcParams::M; m++) {
//				//ippsSum_64fc(sfDataPos, QmcParams::NumP, sfSumPos);
//				(*sfSumPos) = accumulate(sfDataPos, sfDataPos + QmcParams::NumP,
//						complex<double>(0., 0.));
//				//sfSlice += std::norm(*sfSumPos);
//				//(pow((*sfSumPos).re, 2) + pow((*sfSumPos).im, 2));
//				sfDataPos += QmcParams::NumP;
//				sfSumPos++;
////				ippsSum_64fc(sfDataPos, QmcParams::NumP, sfSumPos);
////				sfDataPos += QmcParams::NumP;
////				sfSumPos++;
//			}
//			fftPlan_.execute();
//			//ippsDFTFwd_CToC_64fc(sfSum_, sfSum_, ppDFTSpec_, fftBuffer_);
//			sfPow_[0] += AvNumOfKVev * std::norm(sfSum_[0]);
//			//		* (pow((*sfSum_).re, 2) + pow((*sfSum_).im, 2));
//			sfSumPos = sfSum_.begin();
//			for (int m = 1; m < QmcParams::M / 2 + 1; m++) {
//				double posPS = std::norm(sfSumPos[m]);
//				double negPS = std::norm(sfSumPos[QmcParams::M - m]);
//				sfPow_[m] += 0.5 * AvNumOfKVev * norm * (posPS + negPS);
//			}
//		}
//		*(itObs++) << sfPow_;
//	}
//}
//#if SIM_DIM==2
//void Measure::measureDensityPhonon() {
////	static int i = 0;
////	if (true) {
////		String file_name = QmcParams::outPath + "/WL" + boost::lexical_cast
////				< String > (i) + ".h5";
////		WorldLines::dumpToFile(file_name);
////	}
//	int kmax = int(QmcParams::period);
//	double dp = 1. / QmcParams::period;
//	std::vector<BeadPos> sortTemp(QmcParams::NumP);
//	std::list<BeadPos> notFoundPos;
//	std::list<IntArray> notFilledLatPoints;
//	Array2D<bool> checkFilled(kmax, kmax);
//	Array2D<BeadPos> particlePos(kmax, kmax);
//	for (int m = 0; m < QmcParams::M; m++) {
//// NotFilles
//		std::fill(checkFilled.begin(), checkFilled.end(), false);
//		WL::getWLSlice(m, sortTemp);
//		notFoundPos.clear();
//		std::vector<BeadPos>::iterator itSortTemp = sortTemp.begin();
//		for (; itSortTemp != sortTemp.end(); itSortTemp++) {
//			IntArray indexC;
//			BeadPos &beadPosNew = *itSortTemp;
//			for (int d = 0; d < QmcParams::dim; d++) {
//				indexC[d] = round(beadPosNew[d] * double(kmax)) + kmax / 2;
//				if (indexC[d] == kmax)
//				indexC[d] = 0;
//			}
//			bool isFilled = checkFilled(indexC);
//			if (isFilled) {
//				BeadPos beadPosOld = particlePos(indexC);
//				BeadPos latPos, delta;
//				double distO, distN;
//				std::transform(indexC.begin(), indexC.end(), latPos.begin(),
//						[&kmax](int
//								pos) {return double(pos-kmax/2)/double(kmax);});
//				latPos.subBeadPos(delta, beadPosOld);
//				distO = delta.norm2();
//				latPos.subBeadPos(delta, beadPosNew);
//				distN = delta.norm2();
//				if (distO > distN) {
//					//  Switch
//					particlePos(indexC) = beadPosNew;
//					notFoundPos.push_back(beadPosOld);
//				} else
//				notFoundPos.push_back(beadPosNew);
//			} else {
//				checkFilled(indexC) = true;
//				particlePos(indexC) = beadPosNew;
//			}
//		}
//// Find
//		notFilledLatPoints.clear();
//		for (int kx = 0; kx < kmax; kx++)
//		for (int ky = 0; ky < kmax; ky++) {
//			IntArray index = {kx, ky};
//			if (!checkFilled(index))
//			notFilledLatPoints.push_back(index);
//		}
//// std::cout << particlePos << std::endl;
//		std::list<IntArray>::iterator itNotFilles = notFilledLatPoints.begin();
//		for (; itNotFilles != notFilledLatPoints.end(); itNotFilles++) {
//			IntArray &latIndex = *itNotFilles;
//			double minDist = 10;
//			std::list<BeadPos>::iterator itMinBead, itBead;
//			for (itBead = notFoundPos.begin(); itBead != notFoundPos.end();
//					itBead++) {
//				BeadPos delta, latPos;
//				double distC;
//				BeadPos &currentPos = *itBead;
//				std::transform(latIndex.begin(), latIndex.end(), latPos.begin(),
//						[&kmax](int
//								pos) {return double(pos-kmax/2)/double(kmax);});
//				latPos.subBeadPos(delta, currentPos);
//				distC = delta.norm2();
//				if (distC < minDist) {
//					minDist = distC;
//					itMinBead = itBead;
//				}
//			}
//			particlePos(latIndex) = *itMinBead;
//			notFoundPos.erase(itMinBead);
//		}
//		for (int kx = 0; kx < kmax; kx++)
//		for (int ky = 0; ky < kmax; ky++) {
//			BeadPos delta, latPos;
//			IntArray latIndex = {kx, ky};
//			std::transform(latIndex.begin(), latIndex.end(), latPos.begin(),
//					[&kmax](int
//							pos) {return double(pos-kmax/2)/double(kmax);});
//			particlePos(kx, ky).subBeadPos(delta, latPos);
//			for (int d = 0; d < QmcParams::dim; d++) {
//				Array3D<double> & deltauRef = deltau_[d];
//				deltauRef(kx, ky, m) = delta[d] + 0.5 * delta[(d + 1) % 2];
//			}
//		}
//	}
//	double normalize = 1 / (double(kmax * kmax * QmcParams::M));
//	for (int d = 0; d < QmcParams::dim; d++) {
//		Array3D<std::complex<double>> & deltauqwRef = deltauqw_[d];
////Array3D<double> & deltauPSRef = sfPhononqw_[d];
//		Array3D<double> & deltauRef = deltau_[d];
////	deltauRef.dumpToFile("/home/snirgaz/sim_data/A1.h5");
//		Plans_[d].execute();
////	deltauqwRef.dumpToFile("/home/snirgaz/sim_data/B1.h5");
//		auto itObs = sfPhononqwObs_[d].begin();
//		for (auto it = deltauqwRef.begin(); it != deltauqwRef.end(); it++) {
//			*itObs << normalize * norm(*it);
//			itObs++;
//		}
////		}
////	std::transform(deltauqwRef.begin(), deltauqwRef.end(),
////			deltauPSRef.begin(), deltauPSRef.begin(),
////			[](std::complex<double> x,double y) {return y+norm(x);});
//////		for (int m = 0; m < QmcParams::M / 2 + 1; m++)
////			for (int p1 = 0; p1 < kmax; p1++)
////				for (int p2 = 0; p2 < kmax; p2++)
////					deltauPSRef(p1, p2, m) += norm(deltauqwRef(p1, p2, m));
////		//	std::cout << deltauPSRef;
//	}
//}
//#endif
//void Measure::measureRho() {
//	double arg;
//	BeadPos delta, jump, beforeJump;
//	jump = WL::jumpPos();
//	beforeJump = WL::getDataPos(WL::worm(), WL::tau());
//	jump.subBeadPos(delta, beforeJump);
//// Rho_->at(loc[0] * QmcParams::DenRes + loc[1]) += 1;
//	for (auto it = rhoQ_.begin(); it != rhoQ_.end(); it++) {
//		array<int, SIM_DIM> kVec = it.getIndex();
//		arg = 0;
//		for (int d = 0; d < SIM_DIM; d++)
//			arg += double(kVec[d]) * delta.getD(0);
//		(*it) << arg;
//	}
//}
//void Measure::measureAngle() {
////static int i = 0;
//	vector<double>::iterator wlData = WL::data();
//	double angleqzero = 0, theta;
//	int dimSize = WL::numOfBeads();
//	vector<complex<double>>::iterator sfDataPos, sfSumPos;
//	vector<complex<double>>::iterator anglePos = angle_.begin();
//	BeadPos pos1, pos2, delta;
//	WL::calcNearestNeighbors(nnIndices_, 6);
//	DoublePointerVec wlDim;
//	double norm = double(QmcParams::NumP * QmcParams::NumP);
//	WL::getDimP(wlDim);
////	if (true) {
////		String file_name = QmcParams::outPath;
////		file_name.append("/WL");
////		file_name.append(boost::lexical_cast<String>(i));
////		file_name.append(".h5");
////		WorldLines::dumpToFile(file_name);
////		i++;
////	}
//	complex<double> angle = { 0, 0 };
//	for (int m = 0; m < QmcParams::M; m++) {
//		for (int n1 = 0; n1 < QmcParams::NumP; n1++) {
//			pos1 = WL::getDataPos(n1, m);
//			for (int nn = 0; nn < 6; nn++) {
//				int nnPos = nnIndices_[m][n1][nn];
//				pos2 = WL::getDataPos(nnPos, m);
//				pos1.subBeadPos(delta, pos2);
//				delta.rotate_inv();
//				theta = atan2(delta[1], delta[0]);
//				angle = std::polar((1. / 6.), 6 * theta);
//				//angle.re += (1. / 6.) * cos(6 * theta);
//				//angle.im += (1. / 6.) * sin(6 * theta);
//			}
//			*anglePos = angle;
//			angle = complex<double>(0, 0);
//			anglePos++;
//		}
//	}
//	for (auto it = angleqw_.begin(); it != angleqw_.end(); it++) {
//		array<int, SIM_DIM> kVec = it.getIndex();
//		for (int i = 0; i < dimSize; i++) {
//			angle_[i] *= polar(1.,
//					2 * PI
//							* (double(kVec[0]) * wlDim[0][i]
//									+ double(kVec[1]) * wlDim[1][i]));
//		}
////		expData_[0:dimSize]=2*PI*(
////				double(kVec[0])*wlDim[0][0:dimSize]
////				+double(kVec[1])*wlDim[1][0:dimSize]);
////		ippsCIS_64fc_A50(expData_, sfData_, dimSize);
////		ippsMul_64fc_I(angle_, sfData_, dimSize);
//		sfDataPos = sfData_.begin();
//		sfSumPos = sfSum_.begin();
//		double sfSlice = 0;
//		for (int m = 0; m < QmcParams::M; m++) {
//			(*sfSumPos) = accumulate(sfDataPos, sfDataPos + QmcParams::NumP,
//					complex<double>(0., 0.));
//			//ippsSum_64fc(sfDataPos, QmcParams::NumP, sfSumPos);
//			(*sfSumPos) /= double(QmcParams::NumP);
//			sfDataPos += QmcParams::NumP;
//			sfSumPos++;
//		}
//		if (kVec[0] == 0 && kVec[1] == 0)
//			for (int m = 0; m < QmcParams::M; m++) {
//				angleqzero += 1 / double(QmcParams::M) * std::norm(sfSum_[m]);
//				//* (pow(sfSum_[m].re, 2) + pow(sfSum_[m].im, 2));
//			}
//		fftPlan_.execute();
//		//ippsDFTFwd_CToC_64fc(sfSum_, sfSum_, ppDFTSpec_, fftBuffer_);
//// Power Spectrum
//		transform(sfSum_.begin(), sfSum_.end(), sfPow_.begin(),
//				[](complex<double> sf)->double {return std::norm(sf);});
//		//ippsPowerSpectr_64fc(sfSum_, &(sfPow_.data()[0]), QmcParams::M);
//		*it << sfPow_;
//	}
//	angleZeroQ_ << (angleqzero);
//}
////void Measure::measureJqJq() {
////	int numOfBeads = WL::numOfBeads();
////	int numOfPoints = 2 * WL::numOfBeads();
////	int numOfMeasurePoints = QmcParams::numJqJqPoints * QmcParams::numJqJqPoints
////			* QmcParams::M;
////	DoublePointerVec vPos(QmcParams::dim), avPos(QmcParams::dim), wlPos(
////			QmcParams::dim);
////// Fix Pointers
////	velocity_->getSliceP(0, vPos);
////	averagePos_->getSliceP(0, avPos);
////	WL::getSliceP(0, wlPos);
////// Calculate average Position
////	(averagePos_->data())[0:numOfPoints]=(WL::data())[0:numOfPoints]+0.5*(velocity_->data())[0:numOfPoints];
////	ippsRealToCplx_64f((velocity_->data()), NULL, velocityComplex_,
////			numOfPoints);
////	for (auto it = angleqw_.begin(); it != angleqw_.end(); it++) {
////		boost::array<int, SIM_DIM> kVec = it.getIndex();
////		expData_[0:numOfBeads]=2*PI*(
////				double(kVec[0])*avPos[0][0:numOfBeads]
////				+double(kVec[1])*avPos[1][0:numOfBeads]);
////		// Calc Cis
////		ippsCIS_64fc_A50(expData_, Jq_, numOfBeads);
////		Jq_[0:numOfBeads]=
////		ippsCopy_64fc(Jq_, Jq_ + numOfBeads, numOfBeads);
////		// Calc J(q)
////		ippsMul_64fc_I(velocityComplex_, Jq_, numOfPoints);
////		Ipp64fc *jqPos = Jq_, *jqSumPos = JqSum_;
////		for (int m = 0; m < QmcParams::M; m++) {
////			for (int d = 0; d < QmcParams::dim; d++)
////				ippsSum_64fc(jqPos, QmcParams::NumP, jqSumPos);
////			jqPos += QmcParams::NumP;
////			jqSumPos++;
////			// Fourier Transform
////			for (unsigned int d = 0; d < QmcParams::dim; d++) {
////				ippsDFTFwd_CToC_64fc(JqSum_ + d * QmcParams::M,
////						JqSum_ + d * QmcParams::M, ppDFTSpec_, fftBuffer_);
////				// Power Spectrum
////				ippsPowerSpectr_64fc(JqSum_ + d * QmcParams::M, JqPow_,
////						QmcParams::M);
////				ippsAdd_64f_I(JqPow_, jqjqPos[d], QmcParams::M);
////				jqjqPos[d] += QmcParams::M;
////			}
////		}
////	}
////}
//void Measure::measureWindingNumber() {
////Ipp64u start = ippGetCpuClocks();
//	DoublePointerVec vPos;
//	DoubleVec sum(QmcParams::dim);
//	velocity_.getSliceP(0, vPos);
//	int numOfBeads = WL::numOfBeads();
//// Calc Winding Number
//	double total = 0;
//	for (int d = 0; d < QmcParams::dim; d++) {
//		sum[d] = accumulate(vPos[d], vPos[d] + numOfBeads, 0.);
//		//ippsSum_64f(vPos[d], numOfBeads, &(sum[d]));
//		total += pow(sum[d], 2);
//	}
//	windSqr_ << total;
//	static int i = 0, totalc = 0;
//	if (i < 5 && total > totalc) {
//		String file_name = QmcParams::outPath + "/WL"
//				+ boost::lexical_cast<String>(i) + ".h5";
//		WorldLines::dumpToFile(file_name);
//		i++;
//		totalc++;
//	}
//}
////windDist_[int(round(sum[d])) + length] += 1;
//
//void Measure::measureEnergy() {
//	double total, ek, ev, ef, av, et;
//	av = double(QmcParams::NumP) * double(QmcParams::M);
//	ek = WLEnergy::getEnergyK();
//	ev = WLEnergy::getEnergyV();
//	ef = WLEnergy::getEnergyF();
//	ek = (double(QmcParams::dim) / (2. * QmcParams::epsilon))
//			- (1 / (av * 4. * QmcParams::lambda * pow(QmcParams::epsilon, 2)))
//					* WLEnergy::getEnergyK();
//	ev = (1 / QmcParams::epsilon) * QmcParams::greenFactorV
//			* WLEnergy::getEnergyV() / av;
//	ef = WLEnergy::getEnergyF() * (6. / 9.) * pow(QmcParams::epsilon, 2)
//			* (QmcParams::lambda / av);
//	et = ek + ev + ef;
//	EK_ << ek;
//	EPV_ << ev;
//	EPF_ << ef;
//	E_ << et;
//}
//bool Measure::performMeasure() {
//	OC open_close = WL::getOpenClose();
//	if ((open_close == CLOSED) && (numOfClosed_ < QmcParams::minNumOfMeasure)) {
//		assert(numOfClosed_ < QmcParams::minNumOfMeasure);
//		if (QmcParams::toMeasure.Winding || QmcParams::toMeasure.JqJq)
//			WL::calcVelocity(velocity_);
//		if (QmcParams::toMeasure.Density)
//			this->measureDensity();
//		if (QmcParams::toMeasure.DensitySym)
//			this->measureDensitySym();
//#if SIM_DIM==2
//		if (QmcParams::toMeasure.DensityPhonon)
//		this->measureDensityPhonon();
//#endif
//		if (QmcParams::toMeasure.Winding)
//			this->measureWindingNumber();
////		if (QmcParams::toMeasure.JqJq)
////			this->measureJqJq();
//		if (QmcParams::toMeasure.Energy)
//			this->measureEnergy();
//		if (QmcParams::toMeasure.Angle)
//			this->measureAngle();
//		numOfClosed_++;
//		++(*(measureProgress_.get()));
//		return true;
//	} else if (numOfOpen_ < QmcParams::minNumOfMeasure) {
//		if (QmcParams::toMeasure.Rho)
//			this->measureRho();
//		numOfOpen_++;
//		++(*(measureProgress_.get()));
//		return true;
//	}
//	return false;
//}
//void Measure::dumpMeasurements(String fileName, std::ofstream &outPutFile) {
//	String dsName;
//// Open File
//	hdf5RW hdf5File;
//	hdf5File.initNew(fileName);
//// SF
//	std::cout << "Dumping Measurements" << std::endl;
//	if (QmcParams::toMeasure.Density) {
//		std::cout << "SF" << std::endl;
//		sf_.dumpToFile(hdf5File, "/SF");
//		std::cout << "SF QW" << std::endl;
//		sfqw_.dumpToFile(hdf5File, "/SF_QW");
//	}
//#if SIM_DIM==2
//	if (QmcParams::toMeasure.DensityPhonon) {
//		for (int d = 0; d < QmcParams::dim; d++) {
//			dsName = "SF_PHONON" + boost::lexical_cast < String > (d);
//			sfPhononqwObs_[d].dumpToFile(hdf5File, dsName);
//		}
//	}
//#endif
//
////WindSqr
//	if (QmcParams::toMeasure.Winding) {
//		std::cout << "WIND" << std::endl;
//		windSqr_.dumpToFile(hdf5File, "/WindSqr");
//	}
//// JqJq
//	if (QmcParams::toMeasure.JqJq) {
////		hsize_t dims[4] = { SIM_DIM, QmcParams::numJqJqPoints,
////				QmcParams::numJqJqPoints, QmcParams::M };
////		hsize_t cdims[4] = { 1, 1, 1, QmcParams::M };
////		dsName = "/JqJq";
////		hdf5File.write_double_comp(ds_name, 4, dims, cdims, JqJq_);
//	}
//// Energy
//	if (QmcParams::toMeasure.Energy) {
//		std::cout << "ENERGY" << std::endl;
//		EK_.dumpToFile(hdf5File, "/EK");
//		EPV_.dumpToFile(hdf5File, "/EPV");
//		EPF_.dumpToFile(hdf5File, "/EPF");
//		E_.dumpToFile(hdf5File, "/ENERGY");
//	}
//// rhoQ
//	if (QmcParams::toMeasure.Rho) {
//		std::cout << "RHO" << std::endl;
//		rhoQ_.dumpToFile(hdf5File, "/RHO_Q");
//	}
//	if (QmcParams::toMeasure.Angle) {
//		std::cout << "ANGLE" << std::endl;
//		angleqw_.dumpToFile(hdf5File, "/ANGLEQZERO");
//		angleZeroQ_.dumpToFile(hdf5File, "/ANGLE_QW");
//	}
//	if (QmcParams::toMeasure.DensitySym) {
//		std::cout << "Density Sym" << std::endl;
//		sfqwSym_.dumpToFile(hdf5File, "/DENSITYSYM");
//	}
//	std::cout << "NUM_OF_CLOSED" << std::endl;
//	hdf5File.writeSingle<long int>("/NUM_OF_CLOSED", numOfClosed_);
//	hdf5File.writeSingle<long int>("/NUM_OF_OPEN", numOfOpen_);
//}
//void Measure::readMeasurements(String fileName, std::ofstream &outPutFile) {
//	String ds_name;
//// Open File
//	int numClosed, numOpen;
//	hdf5RW hdf5File;
//	hdf5File.initRead(fileName);
//	hdf5File.read<long int>("/NUM_OF_CLOSED", &numOfClosed_);
//	hdf5File.read<long int>("/NUM_OF_OPEN", &numOfOpen_);
//// SF
//	if (numOfClosed_ > 0) {
//		if (QmcParams::toMeasure.Density) {
//			sf_.readFromFile(hdf5File, "/SF");
//			sfqw_.readFromFile(hdf5File, "/SF_QW");
//		}
//		if (QmcParams::toMeasure.DensitySym) {
//			sfqwSym_.readFromFile(hdf5File, "/DENSITYSYM");
//		}
//#if SIM_DIM==2
//		for (int d = 0; d < QmcParams::dim; d++) {
//			String dsNamePhonon = "SF_PHONON" + boost::lexical_cast < String
//			> (d);
//			sfPhononqwObs_[d].readFromFile(hdf5File, dsNamePhonon);
//		}
//#endif
////WindSqr
//		if (QmcParams::toMeasure.Winding) {
//			windSqr_.readFromFile(hdf5File, "/WindSqr");
//		}
//// JqJq
////		if (QmcParams::toMeasure.JqJq) {
////			JqJq_.
////			hdf5_rw::read_double(fileName, "/JqJq", JqJq_);
////#if SIM_DIM==2
////			int sizesfwq = QmcParams::numJqJqPoints * QmcParams::numJqJqPoints
////					* QmcParams::M * QmcParams::dim;
////#endif
////#if SIM_DIM==3
////			int sizesfwq = QmcParams::numJqJqPoints * QmcParams::numJqJqPoints
////			* QmcParams::numJqJqPoints * QmcParams::M * QmcParams::dim;
////#endif
////JqJq_		[0:sizesfwq]*=double(numOfClosed_);
////	}
//// Energy
//		if (QmcParams::toMeasure.Energy) {
//			EK_.readFromFile(hdf5File, "/EK");
//			EPV_.readFromFile(hdf5File, "/EPV");
//			EPF_.readFromFile(hdf5File, "/EPF");
//			E_.readFromFile(hdf5File, "/ENERGY");
//		}
//		if (QmcParams::toMeasure.Angle) {
//			angleZeroQ_.readFromFile(hdf5File, "/ANGLEQZERO");
//			angleqw_.readFromFile(hdf5File, "/ANGLE_QW");
//		}
//	}
//	if (numOfOpen_ > 0) {
//// rhoQ
//		if (QmcParams::toMeasure.Rho) {
//			rhoQ_.readFromFile(hdf5File, "/RHO_Q");
//		}
//	}
//}
//bool Measure::stopMeasure() {
//	bool flagClosed = ((numOfClosed_ >= QmcParams::minNumOfMeasure)
//			|| (!closedMeasure_));
//	bool flagOpen = ((numOfOpen_ >= QmcParams::minNumOfMeasure)
//			|| (!openMeasure_));
//	return flagClosed && flagOpen;
//}
//int Measure::getNumMeasure() {
//	return std::min(numOfClosed_, numOfOpen_);
//}
//int Measure::getNumClosedMeasure() {
//	return numOfClosed_;
//}
//int Measure::getNumOpenMeasure() {
//	return numOfOpen_;
//}
//void Measure::setNumClosedMeasure(int numMeasure) {
//	numOfClosed_ = numMeasure;
//}
//void Measure::setNumOpenMeasure(int numMeasure) {
//	numOfOpen_ = numMeasure;
//}
//void Measure::initializeCountOutput(std::ofstream &outPutFile) {
//	measureProgress_.reset(
//			new boost::progress_display(QmcParams::minNumOfMeasure * 2,
//					outPutFile));
//	(*(measureProgress_.get())) += numOfOpen_ + numOfClosed_;
//}
//
